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A reaction network is a chemical system involving multiple reac- 
tions and chemical species. Stochastic models of such networks treat 
the system as a continuous time Markov chain on the number of 
molecules of each species with reactions as possible transitions of 
the chain. In many cases of biological interest some of the chemical 
species in the network are present in much greater abundance than 
others and reaction rate constants can vary over several orders of 
magnitude. We consider approaches to approximation of such mod- 
els that take the multiscale nature of the system into account. Our 
primary example is a model of a cell's viral infection for which we ap- 
ply a combination of averaging and law of large number arguments to 
show that the "slow" component of the model can be approximated 
by a deterministic equation and to characterize the asymptotic dis- 
tribution of the "fast" components. The main goal is to illustrate 
techniques that can be used to reduce the dimensionality of much 
more complex models. 

1. Stochastic models for reaction networks. A reaction network is a 
chemical system involving multiple reactions and chemical species. The sim- 
plest stochastic model for a network treats the system as a continuous time 
Markov chain whose state X is a vector giving the number of molecules of 
each species present with each reaction modeled as a possible transition for 
the state. The model for the kth reaction is determined by a vector of in- 
puts z^fc specifying the number of molecules of each chemical species that are 
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consumed in the reaction, a vector of outputs i^'f, specifying the number of 
molecules of each species that are created in the reaction, and a function of 
the state Afc(x) that gives the rate at which the reaction occurs. Specifically, 
if the reaction occurs at time t, the new state becomes 

X{t)=X{t-)+iyi-iyk. 

Let Rkit) denote the number of times that the kth reaction occurs by time t. 
Then the state of the system at time t can be written as 

x{t) = x{o) + Y,Rk{Wk - ^k) = w - '^)m, 

k 

where v' is the matrix with columns given by the z/^, u is the matrix with 
columns given by the z^^, and R{t) is the vector with components Rk{t). 

Rk is a counting process with intensity Xk{X{t)) (called the propensity in 
the chemical literature) and can be written as 

Rk{t)=Yk(^j\k{X{s))ds'^ 

where the are independent unit Poisson processes. Note that writing iJ^ 
in this form makes it clear why is referred to as a rate. 

Defining |i^fc| = X]j ^ik, the stochastic form of the law of mass action says 
that the rate should be given by 

xN/ ^ Ui'^ik^- ( ^ \ AT riii^ifc! ( X 



'iVkfcl-l \iyik ■ ■ ■ I'nik J ' N^^kl yflk-'-Vrnk, 

where is a scaling parameter usually taken to be the volume of the system 
times Avogadro's number, and Kk is a constant specifying the rate of the 
reaction. Note that the rate is proportional to the number of distinct subsets 
of the molecules present that can form the inputs for the reaction. Intuitively, 
this assumption reflects the idea that the system is well stirred, in the sense 
that all molecules are equally likely to be at any location at any time. 

1.1. Law of large numbers and diffusion approximations. If is the vol- 
ume times Avogadro's number and x gives the number of molecules of each 
species present, then c = N~^x gives the concentrations in moles per unit 
volume. With this scaling and a large volume (where large can be pretty 
small since Avogadro's number is 6 x 10^^), 

(1.1) X^{x)^NKu\{c^^'^N\k{c). 

i 

Since the law of large numbers for the Poisson process implies N~^Y{Nu) w 
u, (1.1) implies 

C{t) = N~'X{t) ^ C{0) + Y. f^kli C(s)r {u'f, -Uk)ds, 
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which in the large volume limit gives the classical deterministic law of mass 
action 

k i 

Similarly, since an appropriately renormalized Poisson process can be ap- 
proximated by a standard Brownian motion, that is. 



N 

we can derive a diffusion approximation for the Markov chain by replacing 
Yk{Nu) by y/NWkiu) + Nu, that is, 

+ f'F{C''{s))ds, 
Jo 

where 

F(c) = 5^Afc(c)K-i.fc). 

k 

The diffusion approximation is given by the equation 

C'^(t) = (7^(0) + J2 N-'/^Wk XkiC^'is)) ds^ {u'^ - vk) 

+ f'FiC^{s))ds, 
Jo 

which is distributionally equivalent to the Ito equation 

C^{t) = C^(0) + Y: f ^\k{CHs))dWk{s)K - u,) 
k -"^ 

+ ['F{C^{s))ds 
Jo 

= (7^(0) + yiV"i/2 f\^c^^s))dW 
k 

+ f F{C''{s))ds, 
Jo 
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where a{c) is the matrix with columns y Afc(c)(z^[, — f^). A precise version 
of this approximation is given in [7]. (See also [4], Chapter 10, [5], Chapter 
7, and [12].) 

1.2. Multiscale approximations. Interest in modeling chemical reactions 
within cells has led to renewed interest in stochastic models, since the num- 
ber of molecules involved, at least for some of the species, may be sufficiently 
small that the deterministic model does not provide a good representation 
of the behavior of the system. Modeling is further complicated by the fact 
that some species may be present in much greater abundance than others. In 
addition, the rate constants may vary over several orders of magnitude. 
With these two issues in mind, we consider a different approach to deriving 
a scaling limit approximation of the model. 

will still denote a scaling parameter for the model, but it is no longer 
interpreted in terms of volume or Avogadro's number. In fact, A^"^ plays the 
same role as e in a perturbation analysis of a deterministic model (see, e.g., 
[11]). A may have no physical meaning, but it will have a specific (hopefully 
large) value in any physical or biological setting in which the approximation 
is applied. 

For example, let N be of the order of magnitude of the abundance of the 
most abundant species in the system. For each species i, we then specify 
a parameter < < 1 and normalize the number of molecules by A°', 
defining 

Z,it) = N-''^Xi{t). 

Oi should be selected so that Zj = 0(1), but that still leaves a degree of 
arbitrariness regarding the selection. Note that Oi could be zero, so Zj could 
still be integer-valued. 

We want to express the reaction rates in terms of Z rather than X and 
also to take into account large variation in the reaction rates. Consequently, 
we introduce another set of exponents (3k for the reactions and now assume 
that the reaction rates can be written as N^''Xk{z), where Afc(z) = 0(1) for 
all relevant values of z. The model becomes 

Z,{t) = Zi{0) + ^ A-'n (^J^ A^''Afc(Z(s)) ds^ {u'k - Uk). 

Our goal is to derive simplified models under the assumption that A is 
large, where "large" may be much smaller than Avogadro's number. We 
demonstrate that this process may lead to interesting and reasonable models 
by analyzing a number of examples in the literature. 
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1.3. Outline of the paper. Reaction networks of interest in biology can be 
very high dimensional involving many chemical species and many reactions. 
Consequently, there has been considerable effort to exploit the multiscale 
nature of these systems to derive reduced models. In Section 2 we borrow 
examples from a number of these papers to illustrate how the kind of scaling 
limits we have in mind can be used to provide a rigorous and intuitive 
approach to model reduction. The primary focus of the paper is a model of 
an intracellular viral infection given in [10] and studied further in [6]. We 
analyze this model in Section 3, and give a systematic identification of the 
scaling parameters in Section A. 2. 



2. Examples. 



2.1. Simple crystallization. We consider a model studied by Haseltine 
and Rawlings [6] using the parameters in Table I of their paper. The system 
involves four species and two reactions: 



2A^B, A + C^D. 



The model satisfies 



XA{t) = Xa{0) - 2Yi ^* \kiXa{s){Xa{s) - 1) ds 

-Y2(^j\2XA{s)Xc{s)ds^, 
Xeit) = Xb{0) + Yi ^kiXa{s){Xa{s) - 1) ds^ , 

Xc{t) = Xc{0) - Y2 (^J^ K2Xa{s)Xc{s) ds^ . 

Following Rawhngs and Haseltine, Xa(,0) = 10^, Xb{0) = 0, Xc{0) = 10, and 
Ki = K2 = 10~^. Let = 10^, and take a a = = 1 and ac = 0. Writing 
Ki = K2 = 0.1 X N~^, the normalized system becomes 

Z^{t) = 1 - N"^2Yi (^N j\.05Z^{s){Z^{s) - N~^) ds 

-N-^Y2^l\.lZ^is)Z^is)ds^, 
Z^ it) = A-iyi (^nJ^ 0.05Z^{s){Z^{s) - N~') ds^ , 
Z^ (t) = 10 - y2 (/^* O.lZf (s)Z^ (s) ds^ . 
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Letting — > oo, the simplified system is 



ZA{t) = l- [\.lZA{sfds, 
Jo 

ZBit)= f\mZAisfds, 
Jo 

Zc{t) = 10 - Fa ^* G.IZa{s)Zc{s) ds 



which gives 



ZA{t) ^ 



1 + O.U 

Since Za is deterministic, Zq is a Unear death process with time-varying 
rate A(t) = O.lZ^(t). Consequently, for any t > 0, the distribution of Zc{t) 
is Binomial(10,p(t)) with 

p(t) = exp{-^*0.1Z.(.)<i.} = :^. 

In particular, 

E[Zc{t)] = -^^, Var[Zc(t)] ^ 



1 + O.lt' ' " (l + 0.1t)2' 

compare favorably with the simulation results in Figure 1 of [6]. 

2.2. Enzyme kinetics. Rao and Arkin [9] analyze a model of enzyme 
kinetics, involving an enzyme, substrate, their enzyme-substrate complex 
and a product of this complex 

E + S^ES, ES^P + E. 

The state of this system can be represented by 

ft 



X,{t) = Xs(0) + y_i ^ K_iXes(r) dr 

-Yi(^j\iX,{r)X,{r)dry 
Xes(t) = Xes(O) - y_i K-lXes(r) dr^ 



+ Yi ( /* KiX,{r)X,{r) dr^ - ¥2 (^J^ K2X,,{r) dr 



X,{t) = Xe(0) + y_i ^* K-lXes(r) dr 
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- Yi (^J^ KiX,{r)Xs{r) dr^ + ¥2 ^* AC2^cs(r) dr^ , 



where Xg gives the number of substrate molecules, Xg the number of en- 
zymes, Xes the number of enzyme complexes, and Xp the number of molecules 
of the reaction product. Following Rao and Arkin, take Xs{0) = 100, ^e(O) = 
1000, Ki = K_i = 1, K2 = 0.1. Let N = 1000, and define = N-'^/^X^, 
Z^^ = iV-2/3Xes, Z^ = A^-2/3Xp, = A^'^Xe, and K2 = N-'^/^. Then the 
normalized system becomes 



-Ar-2/3y2^ ArV3^^^(^)rf^ 
(t) = 1 + N-'Y., ^J^ N^r^Zg (r) dr) 



Zp^(t) = Ar-2/3y2^*Ari/3zi^( 



(r) (ir 

Rescaling time by N'^/^ and defining V^^ {t) = Zf{N^/^t), we have 
V^^{t) = 1 + Af-2/3y_^ (^^* Ary^f (r) dr^ 

- Ar-2/3y, (^^* ArVe^(r)K^(r) dr) , 

(t) = -A^-2/3^_i^* Ary,f (r)dr) 
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+ N-^/^Yi (^J^ N^Vj^ {r)Vf (r) dr 

= 1 + iV-iy_i iVKf (r) dr 



- N-^Yi (^J^ N^Vf{r)Vf{r) dr^ 
+ N~^Y2^J^N^/^V,^{r)dr^, 
V^^{t) = N~'/'Y,^l\'/'V^^{r)dr^, 

Vl'it) + K^it) = 1 - N~'/'Y, N'/'K^ir) dr) . 

Note that Vj^ + V;^ < 1 and T^^ + N^^^^Vj^ = 1, so sup^,<i \Vf (r) - 1| ^ 0. 
It follows that for < /3 < 1 and sufficiently large, infj.<t V^^(r) > p and 
V/^it) < Vj^it), where Vj^ is the solution of 

V;^(t) = 1 + iV"2/3y_^ (^J^ NV^!^{r) dr^ - N-^/^Yi ^* N^pVf{r) dr^ . 

The fact that supj.<^ ^^^^p^ ~^ ^ ensures that supg^j,^^Vj^ (r) — > 0, for < 
6 <t, and {V^ , Vj^ + V^) converges to the solution of 

Vp{t)= /Ves(r)dr, 





VUt) = l- tVcs{r)dr, 





that is, Vcs{t) = e~* and Vp{t) = 1 — e~*. On the original time scale (t) w 
1 — e"*/^'', that is, Xp{t) ~ 100(1 — 6"*/"'^'^), which matches well the simulation 
results in the lower plot in Figure 1 of [9]. 

Rao and Arkin also consider Xs{0) = 100, -'^c(O) = 10, ni = = 1, K2 = 
0.1. For this example, let N = 100 and define = N-'^X^, Z^^ = N'^/'^X^^, 
= N-^/'^Xc, and Z^ = Xp, and set K2 = iV*^/^. Then the normalized 
system becomes 

(t) = 1 + N-^Y^i N^I^Zg (r) dr^ - N'^Y^ ^* N^''' Z^ {r)Z^ {r) dr^ 
Zg{t) = -N-^/^Y^^I^jy^/^Z^,{r)dr^+N-^'^Y^I^jy-'/^Z^{^^ 
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(r) dr 



- N-'/^Yi (^Ij N^/^Z^{r)Z,''{r) dr^ + N-'/^Y2 {^f^ Zg (r) dr) , 
Z^{t)=Y^(^f^Zg{r)dry 



Since 



N-^'^Yi(^j^ N''^Z^{r)Z^{r)dr 

< 1 + iV-i/2y_i N^I^Zg (r) dr 



+ 



and Zes{t) < 1, it follows that supj,<^ \Z^ (r) — 1| ^ 0, for each t > 0. Not- 
ing that /q Z^{r) dr = t — Z^{r) dr, we must have /q Z^{r) dr ^ t and 

P 



Z^{t)^Y2{t). 



2.3. Reversible isomerization. Next, we consider a model of reversible 
isomerization studied by Cao, Gillespie and Petzold [3]. The model involves 
three chemical species and two reactions: 



1 < '-'2, '-'2 > '-'S- 



The model is given by 





= Xi{0)-Yi 




= X2{0) + Yi 




-""ill''- 








= Xs{0) + Ys 



The first set of parameter values in (34) of [3] give Xi{0) = 1200, ^2(0) = 
600, Xs{0) = and ki = 1, K2 = 2, K3 = 5 x 10"^ Let N = 1000, and define 
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Zf = N~^Xi,Z^ = iV-^Xa, = X3, and K3 = bN^^I^. Then the normal- 
ized system becomes 

(i) = (0) - N^^Yi ^* NZ^{s) ds^ + N-^Y2 (^J^ 2NZ^{s) ds^ , 
Z^{t) = Z^{Q) + N-^Y^ NZ^{s) ds) 

- N-^Y2 ^* 2NZ^{s) ds^ - N^^Y^ ^* 5N"^/^Z^{s) ds^ , 

Z^{t)=Y,[^j\N-^I^Z^{s)dsy 

Assuming that (Zf (0), (0)) ^ (Zi(0), ^2(0)) (which gives Zi(0) = 1.2 
and ^2(0) = 0.6 for the particular values in [3]), the limiting system is 

Zi(t) = Zi(0)- f Zi{s)ds+ f\z2{s)ds, 
Jo Jo 

Z2(t) = Z2(0)+ f Zi{s)ds- [\z2{s)ds, 
Jo Jo 

Z3(t)=0. 

Consequently, 

Zi(t) + Z2(t) = Zi(0) + Z2(0), 

D{t) = Zi{t) - 2Z2(t) = D{0) - 3 /* D{s) ds, 

Jo 

so 

Ziit) = i(Zi(0) - 2Z2(0))e-3* + |(Zi(0) + ^2(0)), 
Z2(t) = -i(Zi(0) - 2Z2(0))e"3* + i(Zi(0) + Z2(0)) 

and 

^lim Z2(t)) = (|, i)(Zi(0) + Z2(0)). 

Defining = Z^{N'^/^t), the system becomes 

[/f (t) = f/f (0) - iV-iyi (^^* N^/'^^U^is) ds^ + N~^Y2 ^* 2iV^/3[/2^(s) ds^ , 

t^2^(*) = Ui'iO) + N-^Yi (^J^ iV^/^^f (s) ds^ - N-^Y2 ^* 2iV5/3^2^(g) 

-7V-iy3(^*5C/2'^(s)ds), 

[/3^(t) = y3(/^*5C/2'^(s)(is 
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and, hence, 

f/f (t) + Ui'it) = C/f (0) + Ui'iO) - N-'Y-s^l^ 5Ui'is) ds ) . 
Dividing the equation for f/j'^ by N'^/'^ ^ it follows that 

lim (fu^{s)ds- f\ui^{s)ds) =0 

N^oo \Jq Jq J 

and, hence, assuming lim7v_»oo(f^i^(0) + 1^2^ (0)) = C, 

lim / U^{s)ds = \Ct 



N- 

tN 



and C/3 converges to 

U3{t) = Ys 



5Ct\ 



that is, a Poisson process with parameter 5C/3. 

Rescaling time by N^^^, and defining V/^ (t) = N^'^Xi{N^/h), we have 

Vf{t) = ^1^(0) - N-'Yi (^J^ N'^/^Vfis) ds^ + N-'Y2 2N'/^Vi'{s) ds^ , 
Vf{t) = Vi'iO) + N-'Yi ^* N'/^Vf{s) ds^ - N-'Y2 2N'/^Vi'{s) ds^ 

-N~'Y3(^l\NVi'{s)ds^, 
Vi'it) = N-'Ysl^J^ ^NVi'is) ds^ . 

Setting R^{t) = V/^{t) + Vi^it) and assuming lim7v_oo -R^(O) = R{0), we 
have iR^,Vf) (R^Vs) satisfying 

R{t) = R{0) - f i?(s) ds, V3{t) = f i?(s) ds, 

which gives 

R{t) = i?(0)e-(5/3)t^ yg(^) ^ ^(o)(l - e-(5/3)t)^ 

For 6>0, 

hm sup {\Vf{r) - lR{t)\ + \Vf {r) - = 0. 

Note that the simulation results given in Figure 1 of [3] appear to be plots of 
{(X3(tfc)(5, Xj(tfc))}, for some 5>0, rather than of {(ti^,Xi(tk))}, where the 
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{tfc} are the jump times of X^. Consequently, their plots show linear decay 
rather than exponential decay. 

Cao, Gillespie and Petzold also study a second set of parameter values 
(35) in [3] , taking Xi (0) = 2000, X2 (0) = X3 (0) = and ki = 10, k;2 = 4 x 10^ , 
K3 = 2. Letting iV = 10^, we now define = N~^Xi,Z^ = X2, = X3. 
The normalized system becomes 



Z^{t) = 0.2 - N^^Yi ^ WNZf{s) ds^ + N-^Y2 (^J^ 4.NZi^{s) ds^ 

Z2{t) = Yi ^* lOiVZf (s) ds^ - Y2 mZ^{s) ds^ - Ys IZ^^is) ds^ , 

Zi'{t) = Y3(^j\zi'{s)ds^. 

Dividing the equation for Z2 by N, we see that 

sup ^ (yi £ WNZ^is) ds^ -Y2 ANZi' {s) ds^ ) ^ 0, 
for each t>0, and hence, 

sup|Zf (r) -0.2| ^0, 

r<t 



sup 

r<t 



Z^{s)ds-^r 



0, 



and Z^ converges to Ys{t). 

Rescaling time by A^, and defining Vj^{t) = N~^Xi{Nt), for z = 1,3 and 
V2^it) = X2{Nt), we have 

V^^{t) = 0.2 - N~^Yi ^* lOAfVi^(s) ds^ + N~^Y2 (^j\N^Vi^ {s) ds 
Vf{t) = Yi ^* lON^Vfis) ds^ - Y2 AN^V2^{s) ds 

-Ysl^j\NV2''is)ds^, 
Vf{t) = N-'Ys ^* 2iVy2^(s) ds] . 



Let V^{t) = Vf{t) + N-'^Vi^{t), and since earlier results imply {t) 

Z^{Nt) PS i, we have sup,,<t \Vl^ {r) - i\^(r)| ^ 0. 

Then 

Vl^{t) = 0.2 - iV-iyg 2NVi^{s) ds 



and from the equation for V2 

tN 
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sup 

r<t 



lOVf {s)ds- / {A + 2N~')V^\s) 
Jo 



Consequently, {Vi , V3 ) converge to the solution of 



=0.2 - f\vi{s)ds, 
Jo 



V3{t)= / 5Viis)ds, 
Jo 



giving 



= 0.2e 



-5t 



FsW =0.2(1 -e 



-5t\ 



To better understand the behavior of , for a bounded function /(f2), 
define 

A^f{vuV2) = Wvi{f{v2 + 1) - fiv2)) + (4 + 2N~')v2{f{v2 - 1) - fiv2)) 
and note that 



f{Vi'{t))-N-' I A'y{Vfis),V^''is))ds 
Jo 

is a martingale. Dividing by N'^, it follows that 



rN I 



rN I 



A^f{V{^{s),Vi'{s))ds^Q 



rN 



N I 



and that for each t, V2 {t) converges in distribution to a random variable 
V2{t) satisfying 



where 



E[Af{V^{t),V2m = ^, 



Af{vi,V2) = mvi{f{v2 + 1) - f{v2)) + Av2{f{v2 " 1) " f {V2)) 



(see, e.g., [8]). 

For each vi, Ay^f{v2) = Af{vi,V2) is the generator of an infinite-server 
queueing model with arrival rate lOvi and service rate 4. It follows that 
V2{t) has a Poisson distribution with parameter 2.514 (t). Note that 
does not converge in a functional sense. In particular, for < ti < t2 < 
••• < tm, {y^i^ {ti),V^ {t2), ■ ■ ■ {tm)) converges in distribution and the 
components of the limit {V2{ti),V2{t2), ■ ■ ■ , V2{tm)) are independent Poisson 
random variables. 
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3. Intracellular viral kinetics. Next we consider a model of an intra- 
cellular viral infection given in [10] and studied further in [6]. We follow 
the presentation (in particular, the indexing) in [6]. The model includes 
three time-varying species, the viral template, the viral genome and the vi- 
ral structural protein. We denote these as species 1, 2 and 3, respectively, 
and let Xi(t) denote the number of molecules of species i in the system at 
time t. The model involves six reactions, designated (28a)-(28f) in [6]: 

(a) T + stuff^T + G, 

(b) G^T, 

(c) r + stuff^r + 5, 

(d) T^0, 

(e) 5^0, 

(f) G + S^V, 

where "stuff" refers to nucleotides and amino acids that are assumed avail- 
able at constant concentrations. The basic model satisfies 



Xi(t)=Xi(0) + n^ K2X2{s)ds^ -Yd(^J^ K^Xi{s)ds 
X2{t) = X2(0) + Ya ^* KiXi{s) ds^ - n ^* K2X2{S) ds 

-Yf(^j\^X2{s)X^{s)ds^ 
Xs{t) = XsiO) + Y,(^J^ K3X,{s) ds^ - Ye ^* ^5X3(3) ds 
-Yf(^j\^X2{s)X3{s)ds 



Following Haseltine and Rawlings, Xi(0) = 1,X2(0) = X3(0) = 0, while 
the reaction constants from Table III of [6] are as given below. Let N = 1000, 
corresponding to the order of magnitude of the largest reaction constant. 
Then the rate constants can be expressed as follows: 





1 


1 




0.025 


2.5iV-2/3 




1000 


N 


K4 


0.25 


0.25 


K5 


2 


2 


Ke 


7.5 X 10-^ 


0.75A^-5/3 



Note that, for simplicity, we have replaced K5 = 1.9985 by K5 = 2. 

We have identified N with the largest rate constant, but it is also the 

7N 

'2 



order of magnitude of the most abundent species. Writing = Xi, 
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jY-2/3j^2 and = N^^X^, the normalized system with the scaled rate 
constants becomes 



Z^{t) = (0) + Yb (^^ 2.5Z^{s) ds^ - Yd ^ 0.25Zf (s) ds 
Z^{t) = Z^{<d) + 7V-2/3y^ z^{s) ds^ - iV-2/3y^ 2.5^2^(5) 

- 7V-2/3y^ 0.75Z2^(s)Z3^(s) ds) , 
Z^{t) = (0) + A^-i^c (^f^ NZ^{s) ds^ - iV-i 



Yj l\NZ^{s)ds 



We also write for the system with rate constants expressed in terms 
of N , including the superscript only to emphasize the dependence of the 
model on the scaling parameter. 

There is substantial probability that the infection dies out quickly, but if 
Z2 reaches any significant level, the chance becomes negligible. 

To be precise, let 

K^{t)=Ya{^j\^{s)dsy 
and for A: = 1, 2, ... , define 

(3.1) (i^ = \nl{t>Q:K^{t)>k]. 

We have the following results. 

Lemma 3.1 (Probability of infection). For [3^ defined in (3.1), 



lim lim < 00} = 0.75. 



In particular, there exist kjy — > 00 such that limjy^oo P{PkM ^ ~ 0.75. 

Remark 3.2. The theorem essentially gives the probability that a single 
virus successfully infects the cell. The argument, which essentially compares 
the initial stages of the infection to a branching process, is a standard tool 
in the analysis of epidemic models. See, for example, [2]. 

Proof of Lemma 3.1. To understand the initial behavior of the sys- 
tem, consider 

(3.2) X^{t)=Yb(^j\.5N-^/^X^{s)ds^ - Yd(^j\.25X(^ (s) ds^ , 
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(3.3) 



Xi'{t) = l + Ya(^I^X^is)ds^-Yb(^j\.5N-^/^Xi'is)ds 



(3.4) 



Z^{t) = N-'Y, (^Ij NX^{s) ds^ - N'^Y, 

- N-% 0.757V-2/3X2^(s)Z3^(s) ds^ . 



2NZ.^(s)ds 



If the virus production term (reaction /) is dropped from the equation 



for X2 , the resulting system, 



X^{t)=Yk(^j\.5N~^/^X^{s)ds^-Yd(^j\.25X^{s) 



ds), 



Xi'{t) = l + Yj I X 



:f (s) ds^ - Yb ^* 2.5iV-2/3^Af ^ 



determines a continuous time, two-type branching process. It is easy to check 
that the process is supercritical. The "lifetime" of each Type 1 molecule is 
exponentially distributed with parameter 0.25 and the number of Type 
2 molecules created by the zth Type 1 molecule during its lifetime has a 
shifted geometric distribution with expectation 4, that is, 



;>oo 

P{i = k]= 0. 
Jo 



25g-o.2Mg-t 



dt 



1 /4 



A: = 0,1,.... 



k\'' 5V5. 

Each Type 2 molecule is eventually converted to Type 1. Starting with a 
single Type 1 or Type 2 molecule, the probability of extinction is simply the 
probability that 

n 

5„ = l + ^(e^-l) 

i=l 

hits zero for some n > 0, an event with probability 0.25. 

To complete the proof of the lemma, one only needs to check that 

Yfi^j^' ^.7^N-^'^X^{s)Z^{s)ds^^Q 
in probability for each k. But 



E 



Yflyj^ ' ^.lbN-^/^X^{s)Z^{s)ds 

<{).7bN-'^/'\k + l)E / Z^{s)ds 
Jo 
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s) ds 



<0.375N''^/''^{k + l)E I " 
<0.375Af~2/^(fc + l)fc 

where the second inequahty follows from equation (3.4) and the last inequal- 
ity follows from the fact that 

r rPi 1 r / /-/^f \i 



E 



We now want to describe the behavior of the process once the infection 
is established. Since we have scaled X2 by N~'^/'^, X2 must reach a level 
that is 0(iV^/^) to have nontrivial behavior. 

As in the proof of Lemma 3.1, if the virus production term is dropped from 
the equation for X2 ^ the expectation of the resulting two- type branching 
process (Xf satisfies m{t) = Q^m{t), where 

-0.25 

1 -2.5iV~2/3 



Q 



N 



Following the classical analysis for branching processes (see Section V.7 
of [1]), the largest eigenvalue for Q^is 



A 



N 



(0.25 + 2.5iV~2/3) + ./(0.25 + 2.5iV-2/3)2 + 7.5iV-2/3 



and the total "population" should grow like e^^*. There exists > sat- 
isfying 

(1 - 0.25p^) = A^p^, 2.5(p^ - l)7V-2/3 = A^, 

that is, {p^ , 1) is the corresponding left eigenvector. It follows that p^ — > 4 
and iV2/3^7V ^7 5^ 

Let R^{t) = p^X^{t)+X^{t), and define 

rf =inf{t:i?^(t)>iV2/3^}. 

We are really interested in the first time X2 reaches N'^^^e, but defining 
in terms of rather than X2 simplifies the proof of the next theorem. 

Theorem 3.3 (Time until establishment). For as in Lemma 3.1 and 
each < e < 2 and 6 > 0, 



lim P 

N-*oo 



iV2/3 log AT 

In particular, liniAr^oo P{t^ < 00} 



4 

45 

0.75. 



>6 



0. 
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Proof. In the calculations that follow, recall that the law of large num- 
bers implies that, for a unit Poisson process, 

Yiu) 



(3.5) 



lim sup 



1 



a.s. 



In addition, note that, without loss of generality, we can assume that /cjv/ 
logiV^O. 



Define 



K^it) = YJ I Xi\s)ds]- I Xi\s)ds, 



Kit) = Y, ^* 2.5N-^/^X^is) ds^ , 

k^{t) = Y^,(^j\.m-^I^X^{s)ds^ - j\.m-^I^X^{s)ds, 

and similarly for , and so on. Since /o"'^ X^{s)ds is the /t^th jump 
time of Ya, 



^ '^X^{s)ds-1 

KN Jo 



0, 



and it follows from (3.5) that 



<oo} ,SUp 



—' k IV 



Kit) 



IoKis)ds 



<oo} sup 



Kit) 



N 



j22.5N-y^X^{s)ds 



and 



l{/3^ <oo} sup 



Kit) 



Kit) 



0.25 



0. 



With reference to (3.4), ETf (i) < K^{t) and 



\ds 



fKit) \ 

[ - 1 VO 



lim 1 



{/3f <oo} sup -fTj^TK 
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Consequently, 



lim l^i^N sup , r^N(+\ - I V 

f^kfj e 



< lim 1 r^AT sup ttJt^ — - - 0.375e V 



< <oo} sup ^^;.V - 0-375e 1 V 

= 0. 

Since K^{t) < K^{t) < 1 + K^{t) and 



R^' it) = l + it) - p^'K^^ it) + (p^^ - l)Ki' it) - Kf it) 

■N 



>1 + K^{t)-K^{t)-Kf{t), 



we have 



(3.6) > 



R^{t) 



> 1 



and 



(3.7) 



1 + K^it) 
K^{t) + Kf{t) 



( K^" [t) + KV it) \ 
lim IrgN inf 1 / \^/vA - 0-75(1 - 0.5e) A 



0. 

In other words, on the event {j3^ < co}, 

K^{t)+Kf{t) 



inf 1 , . 



is asymptotically bounded below by 0.75(1 — 0.5e) > 0. 
Writing p instead of p^ , and (5 instead of , 

logi?^(i) =logi?^(/?fj + /Sog^^^^^dK,^(.) 
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(1 - 0.25p)Xf (a) + - l)2.5Af^^/3j^2^(g) ^ 
_ 0.75iV"^/3X2^(g)Z3^(s) 
The second term on the right-hand-side is bounded by a constant times 



< sup 



and similarly for the third, fourth and fifth terms. By (3.6) and (3.8), 

<oo} sup 

is stochastically bounded and similarly with ^ and . The sixth 
term is a martingale with quadratic variation 

N , 



RN^s-Y 
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and stochastic boundedness of the sequence of martingales follows from the 
stochastic boundedness of the quadratic variation and the boundedness of 
the jumps of the martingale. (See Lemma A.l.) 
The last term satisfies 

<0.75Af-2/3 f'z^{s)ds. 
Jo 



By the equations for , 
E 



^3 



Z^{s)ds 



<0.5E 



X^{s)ds 



and adding the equations for and X2 and taking expectations, 



0.75^ 



+ 1 



= ii;[Xf(rf)+X,^(rf)+Kf(rf)] 

< E[X^iT^) + X2^(rf )] + 0.75eE 

< E[X^{t^) + X^{t^)] + 0.375e£; 



Zi'{s)ds 



X^{s)ds 



and hence, 

0.75(1 - 0.5e)E 

(3.8) 
and 



X^{s) ds 



<i?[Xf(rf)+X,^(r^ 
< iV^/^e + p 



£;[C/^(rf )] < 

Consequently, 

logi?^(Tf) 

= logi2^(/3f^)+0(l) 



2-£ 



(1 - 0.25p)Xf (g) + {p- l)2.hN-'^/^X^{s) 



ds 



:logii^(/?,^J + 0(l) + A^(rf-4^^; 
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Note that 

limsup/^<limsup-— / (Xf (.) + (s)) 
N^oo UnN^i-^ n^oo kmN'^i-^ Jo 

1 



< lim sup irir-K^iPkj^ j 

1 



< 

- 2.5' 



Then assuming kj\f/ log ^ 0, 



N _ oN N 
PkN _ 



A^2/3iogAr Ar2/3iogiv' 
and since {PkJ < (1 + + ^Jv) and A^^/s^JV ^ 7 5^ ^^^^^^ 



rf logi2^(rf) 2 



A^2/3iogiV logA^ 3' 
giving the result. □ 

The computations in the proof of Theorem 3.3 also give the following 
lemma. 

Lemma 3.4. Let < ei < £2. Then N~^/^{t^ - t^) is stochastically 
bounded, and for 5 > 0, 

hm Pi — log — — 6 



N-*oo { 15 £1 

<— log-+/ ^ ^ ^ ds + ,5 = 1. 

15 £1 Jrj^ J 

Remark 3.5. In fact, we will see that N^'^^'^{t^ — rj^) converges to a 
constant. 

Proof of Lemma 3.4. Since for < t < rj^, i?Ar(t) = ©(iV^/^), the 
integral expression for logi?Ar(t) implies 

The lemma follows from the fact that the last term is nonnegative and 
stochastically bounded. □ 
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On t G [0,r^], Xi^{t) is dominated by the linear death process with im- 
migration satisfying 



Xi{t) =Yb{2.5et) -Yd(^j\.25Xi{s)ds 



that is, an infinite server queue with Poisson arrivals of rate 2.5e and expo- 
nential service times with rate 0.25. For /3 > 0, let 7^ = mi{t : Xi{t) > N^^}. 
By Dynkin's formula, for each i > 0, 







7^ At 

' Cf{Xi{s))ds 



E[f{Xi{tAj'^))]=f{0) + E 
where 

Cf{k) = 2.be{f{k + 1) - f{k)) + 0.25kifik - 1) - f{k)). 
For < (5 < 1, let f{k) = {klf . Then cq = sup^ Cf{k) < 00, and 



It follows that, for any a > 0, 

hm P{7^<iV"} = 0. 

Consequently, taking < /3 < 2/3 and a > 2/3, we see that N~^/^X{^{t^) 
and, hence, N-'^/^X^{t^) e. 

It is clear that the time scale for Z2 is much slower than the time scale for 
Zi and Z'^ . The fast time scale of {Zi^,Z^} enables us to "average out" 
their contributions to the evolution of the second component after time . 
We define Vf{t) = Zi{T^ + N^/H). 

Theorem 3.6 (Averaging and deterministic approximation). 

(a) Conditioning on < oo, for each 6 > and t>0, 

lim p( sup 1^2^(5) - V2(s)| > = 0, 
where V2 is the solution of 

(3.9) V2{t)=e+ f\.5V2{s)ds- l\.7W2{sfds. 

Jo Jo 

(b) Conditioning on < 00, for each t > 0, {Vf {t),Vf {t)) converges 
in distribution to a pair {Vi{t),V3{t)) with joint distribution fij^ satisfying 

2.5V2{t){g{z + l,y)-g{z,y)) 



(3.10) 



+ 0.25z{giz - 1, y) - g{z, y)) + {z - 2y)^{z, y) 

oy 



fi',-\dz,dy) = 0. 



24 BALL, KURTZ, POPOVIC AND REMPALA 

In particular, Vi{t) has a Poisson distribution with parameter WV2{t), so 
E[Vi{t)]=YaT{Vi{t)) = 10V2{t), 
E[V3{t)] = 5V2{t), Vav{Vs{t)) = f 

and 

Cov(yi(t),F3(t)) = f^2(t). 

Remark 3.7. (a) Of course, the equation in Part (a) is just the classical 
logistic equation. 

(b) For times ti < ta, [Vi^ {ti),Vf (ti)) and (V^i^(t2), V^a^fe)) are asymp- 
totically independent. 

Proof of Theorem 3.6. On the event <oo, 

Vl'it) = Z^{t^) + ^* 2.5iV2/3yAf _ ye (^J^ 0.25N^/^Vf {s) ds 

Vi^it) = Z^{r^) + N-'/'Y^^J^ N'/'vfis) ds 

- iV-2/3y; (^Ar2/3 o.rbVi" {s^i" (s) ds^ 
Vfit) = Z^{r^) + N-'Y![^f^ N^I^Vfis) ds 

- iv-iy/ ^* 2N^I^Vi^{s) ds^ 

- N-^Yf ( f Q.l^N^'^Vi^ {s)Vi^ {s) ds ] , 



where Y^, Y^ , and so on, are unit Poisson processes obtained from Ya, Yf,, 
and so on. by taking increments over the appropriate intervals. For example, 

Y^^u) = Yb (^j^' 2.5Z^(s) ds + v^- n (^j^' 2.5Zf (s) d^ . 

By the martingale properties of the Poisson processes, 

E[p^N-^/'^Vf{t) + Vi'{t)] 

= E[p^N~^'^Z^{t^) + Z^{t^)] 

+ r ^[(1 - 0.25p^)yi^(s) + 2.5(p^ - l)Vi^{s)] ds 
Jo 
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+ A^iV2/3 /* i?[p^iV"2/3y^^(s) + vi'is)] ds, 
Jo 

so by Gronwall's lemma, 

and the equation for implies 

JO 



and similarly for Vi^ 



rN 

3 ' 



E[Vi'{t)]<E[Zi'{T^)]e~'''''''+ t N^l^e-^''"'^'-'^E[Vf{s)]ds 

Jo 

The law of large numbers for the Poisson process implies that is 
asymptotic to 

Vi'it) = Z^{t^) + [\vf{s) - 2.5Vf{s) - 0.75Vi'{s)Vi'{s))ds, 
Jo 

and since Z2 {t^) £, as in the proof of Lemma A. 4, we can verify the 
relative compactness of {^^} if we can verify the stochastic boundedness of 

f\v!'is) - 2.5Vfis) - 0.75Vfis)Vi'is)f ds, 
Jo 

which in turn will follow from the stochastic boundedness of 

(3.11) fvfis)'^ds, fvi'{s)^ds, fvi'{s)^ds, 

Jo Jo Jo 

for appropriate k. 

Let 7^ = inf{t : Z2 {t) > 4}, and in the equation for Z^ , replace lb(/o 2.5 x 

Z^{s) ds) by YbC/o^^'^ 2.5Z^{s) ds). Note that 7^ > rf , and if we can verify 
the stochastic boundedness of (3.11) for the modified system and show that 
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7^ oo, we will have the stochastic boundedness for the original system. 
Note that 

(3.12) Vf {sf ds = N-^'^ Z^{ufdu. 

Jo Jri^ 

Taking eo = | and to < ^ log 2, it follows from Lemma 3.4 that the sequence 
in (3.12) is stochastically bounded for each t if and only if 

is stochastically bounded for each ti. By Lemma A. 3, 
(3.13) 



'^0 



< {Z^{T,1f V l)e-0-24^'''^ + C7fcll(l - e-0-25A^'''^), 

where the 11 comes from the fact that 2.5Z2 {u)1^^^^n^ + 1 < 11. Since 

{t^^ ) < + 1, the first term goes to zero, and the stochastic bound- 

edness follows. Stochastic boundedness for 



rzt'{r,^^+Ny-'s)'ds, 



i = 2,3, follows by a similar argument, using the fact that (rj^ ) < Eq + 
iV-2/3 and, by (3.8), 



JV 

E[Zi^{r^J]<E r 



< 



0.75(1 -0.5eo)' 

and applying (3.13) to bound the second term on the right-hand-side of 
(A.5). 

As — > 00, dividing the equations for and by N"^'^ shows that 



/* Vf{s) ds - 10 /* Vf{s) ds 0, 
Jo Jo 

t V/(s) - 5 /* Vi^{s) ds 0. 
Jo Jo 



The assertion for V^^ and the fact that V-]^ is asymptotically regular (e.g., 
one can prove that lim^^olinisup^^^oo E[svi\)i<^j' sup5</j \ V2^ {t + s) — V2' {t)\] 
0) implies 



Vi^{s)Vf{s)ds-5 Vi^isy'ds^O. 
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It follows that converges to the solution of (3.9). It should now be clear 
why we shifted the initial time to . 

and fluctuate rapidly and locally in time. behaves like a 
simple birth and death process with entering as a parameter, and 
can be approximated by an ordinary differential equation driven by , 
that is, 

(a + iV-2/3r) « Vf{a)e-^^ + f e-^^'-'WI" {a + N-^'h) ds. 

Jo 

To be specific, let 

B^g{z,y)=2.5Vi'{s){g{z + l,y)-g{z,y)) 
+ 0.25z{g{z-l,y)-g{z,y)) 
+ zN{g{z,y + N~^)-g{z,y)) 
+ 2yN{giz,y-N-')-g{z,y)). 

Then 

M^{t)^g{Vf{t),Vi'{t))-g{Z,{r^),Zs{T^)) 

-N^l^ fB^giVf{s),Vi'{s))ds 
Jo 

is a martingale, and defining an occupation measure by 

r^(Cxi?x [0,t])= f\c{Vl'{s))lD{Vi'{s))ds, 
Jo 

(3.14) (t) =5(Fi^(t),y3^(t)) -5(^i(rf ),Z3(rf )) 

- iV2/3 f B^g{z,y)T^{dz x dy x ds). 

JZ+xR+x[0,t] 

Let Cm = Cm^'^ X IK"*") denote the space of measures v on Z"^ x M+ x 
[0,00) such that z^(Z+ x M"*" x [0,t]) =t, topologized so that convergence is 
weak convergence on x x [0,i] for each t > 0. It is easy to verify that 
the sequence {¥21^^^ is relatively compact in D]g+([0,oo)) x /Z^, where 
-Dig+ ([0, 00)) is the space of cadlag R+-valued functions. Let (V2,r) be a 
limit point of the sequence. 

Define 

Bvd{z, y) = 2.hv{g{z + 1, y) - g{z, y)) + 0.25z(c/(z - 1, y) - g{z, y)) 

+ {z-2y)^{z,y). 
dy 

Noting that 

lim {Bfg{z,y)-ByN,^.g{z,y)) = Q 
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and dividing (3.14) by A^^/^ and letting N oo, we have 
(3.15) / By,,)g{z, y)T{dz x dy x ds) = 0. 

Differentiating (3.15) gives (3.10) at least for almost every t. (See [8] for 
more details.) 

From (3.10), we can easily obtain all the moments of the limiting joint 
distribution. Let {Zs,Ys) be a random vector with the law Note that 
Zs is just a Poisson random variable with expectation 10V2(s). Hence, the 
marginal moments of Zs are immediate. With g{z,y) = y in (3.10), we get 

j[z-2y]^if{dz,dy)=Q 

and, hence, 

E[Ys] = \E[Zs\=W2{s). 

With g{z, y) = zy, we get 

J [z{z - 2y) + lyV2{s) - lyz]f,l^{dz,dy) 



0, 



E[ZsYs] = ^E[Zl] + fV2{s)E[Ys\ = fV2{s) + hm{sf . 
With g{z,y)=y^, 

J[2y{z-2y)]fif{dz,dy)=0 

and 

E[ys] = lE[ZsYs]. 
In general, taking g{z,y) = z'^y'^, one gets the recursive equation 



(3.16) 



n-l 

mE[Zr'Yr'] - 2mE[Z^Yn + IV2{s) ^ ( ^ ) E[zr'Yn 



71-1 ^ V 



from which one can iteratively compute all the joint moments of {Zs,Ys). 
Note that, in order to compute E[Z'^YJ^], one first has to compute 

E[Z^+'^], ElZ^+'^-^Ys], EiZ^+^'-^Ys^ . . . , EiZ^+^Ys"'-^ 

as weh as ah of {E[ZiY^^] :0 < i + j < n + m}. □ 
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APPENDIX 

A.l. Estimates. 

Lemma A.l. Suppose M is a martingale with quadratic variation [M] 
and Z = supt\M{t) - M{t-)\. Then 

p|sup|M(s)| > cl < 1±^[^ + p{[M]t > d}. 

Proof. Let r = inf{s : [M]s > d} and note that [M]rAt <d + Z^. Then 
p(sup|M(s)| >c| <p(sup|M(r As)| >c|+P{r<t}, 
and the result follows by Doob's inequality. □ 
We need the following inequality. 

Lemma A. 2. Let z(t) > 0, t>0, F : [0, oo) — > [0, oo) be nondecreasing, 
and b> 0. Suppose that, for <r <t 

(A.l) z{t)-z{r)< f\F{z{s)) -bz{s))ds. 

Jr 

If z(0) < z*{0) and z* satisfies 



z*{t) = z*iO)+ / iF{z*{s))-bz*{s))ds, 







then z{t) < z*(t), t>0. If, in addition, z*{0) > 1, c > 6, and for some k> 1, 
F{u) < cu^^-^^/^, u>l, then 



z{t) < ('^*(0)i/*^e-(^/'=)* + ^(1 - e-^"/''^') 

(A.2) 



k 



<,*(0)e-(VfcH + i_(i 



e 



Proof. Let zq = z, and let zi satisfy 

(A.3) = z(0) + f\F{z{s)) -bz{s))ds - f\{zi{s) - z{s)) V Ods. 
Jo Jo 

Let T = {r:zi{r) > z{r)}. For r G P, (A.l) and t>r, 
(A.4) z{t)<zi{r)+ l\F{z{s))-bz{s))ds. 
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Let t* = sup{s <t:zi{s) > z{s)}. The continuity of zi and (A. 4) imply that 
z{t*) < zi{t*). U t>t*, then z(s) > zi{s), for t* < s <t, and the last term 
in (A. 3) is zero, so 

z{t) < zi{t*) + J\f{z{s)) - bz{s)) ds = Zl{t). 

Consequently, zi(t) > zo{t), for all t > 0, and 

zi{t) = z{0) + f\F{zo{s)) - bzi{s)) ds. 
Jo 

For m > 0, define Zm+i recursively by requiring 

Zm+l{t) = z{0) + [ {F{Zm{s)) -bZm+l{s))ds. 

Jo 

Then Zm+i > Zm and Zm converges to the solution of 

z*{t) = z{0) + f\F{z*{s)) - bz*{s)) ds. 
Jo 

It follows that z* > z. 

Suppose F{u) < cu^^"^)/*^, u>l,c>b, z*{0) > 1, and z**{t) satisfies 

z**{t) = z*{0) + [\cz**{s)^^-^^'^ - bz**{s))ds. 
Jo 

Note that z**{t) > 1, t > 0, so F{z**{t)) < cz**{t)^^-^^/^ , t > 0, and z*{t) < 
z**{t),t>0. Setting z**{t)=y{t)'', 

y{t) = z*iO)'/^ + ^t- f\y{s)ds 

and 

y{t) = z*(0)Vfce-(^/'=)* + -(1 - □ 



Lemma A. 3. Let a,f3>0. Suppose Z >0 is adapted to a filtration {J-t} 
and 

Z{t) = N-f{Ko + Ki{t) - K2{t) - Ks{t)), 

where Kq is an integer-valued random variable and Ki, K2 and are 
counting processes with no simultaneous jumps and {Tt} -intensities N°^U{t), 
fiN'^Z{t) and AAr(t), respectively. Then, there exists Ck> such that 

E\Z{t)''\Z{0)] < (Z(O)'^ V l)e-^^""''* 

(A.5) 

+ CksnpE[{U{s) + 1)'](1 - e-'^^"""*). 

s<t 
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Proof. Let y be a unit Poisson process that is independent of the Ki, 
and let Z satisfy 

Z{t) = + i^i(t) - K2{t) - Y (^J^ iim{Z{s) - Z{s))ds^ ) . 

Then Z{t) < Z{t) and 

K2{t) = K2{t) + Y(^j\N''{Z{s) - Z{s))d^ 

has intensity fj,N"'Z(t). We have 

Z{t)^ = Z(0)^ + l\{Z{s-) + N-^f - Z{s- f) dKiis) 
Jo 

+ [\{Z{s-)-N-Pf -Z{s-f)dK2{s) 
Jo 

and 

E[z{t)''\zm 

Jo \ ) 

+ ;u(-l)^'-'Z(s)'+iiV-^('="')|Z(0)] ds 
= Z(0)^ + Ar«-/3 f J2E[HN^,^i{U{s),Z{s)\Z{0))]ds 

^0 1=0 

- f' fiN"~'^kE[Z{s)''\Z{0)]ds, 
Jo 

where 

HN,k,i{U{s),Z{s)) 



Applying the Holder inequality, there exists a^,/ > not depending on N 
such that 

E[HN,k,i{U{s),Z{s))\Z{0)] 

< ak,iE[{U{s) + l)''\Zm'/''E[Z{s)''\Z{0)f/\ 

It follows that 

E[z{t)'\zm 
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< Z(O)'^ + iV"-^ r E[{U{s) + lf\Z{0)f^G{E[Z{s)''\Z{0)])ds 
Jo 

' fiN"-f^kE[Z{s)''\Z{0)]ds, 



where G{u) is a polynomial of degree A; — 1 in u^^^ . It follows that, for u > 1, 
there is a constant such that G{u) < afcU^^"-*^)/*"'. Applying (A. 2) with 
b = iiN"~^k and appropriate choice of c gives (A. 5). □ 

Lemma A. 4. Suppose that /? > and that Z^{t) = N~f^K^{t), where 
is a counting process with intensity N^Xi\f{t) . Suppose that {Jq AAr(s)^ ds} 
is stochastically bounded for each t> 0. Then {Z'^} is relatively compact as a 
sequence of processes in the sense of convergence in distribution in D^[0, oo), 
and every limit point has continuous sample paths. 

Proof. Since can be represented as Y{N^ Jq Xiy{s) ds) for a unit 
Poisson process Y, by the law of large numbers, it is enough to verify the 
relative compactness of the sequence A^(t) = /g AAr(s) ds. But for t <t + h < 



\A^{t + h)-A^{t)\<Vh]Jj^XN{syds, 

which gives a uniform equicontinuity condition implying the relative com- 
pactness of {A^} (see, e.g.. Theorem 3.7.2 of [4]). □ 

A. 2. Determining the scaling exponents. The scalings employed for the 
examples in Sections 2 and 3 were determined in part by examining the 
published simulations. In particular, these simulations suggested the rela- 
tionships among the a^. That approach to the choice of the scalings, how- 
ever, is still more art than science and leaves open the question of whether 
slightly different, but equally reasonable scalings would produce significantly 
different limiting approximations. In this section we reconsider the model of 
Section 3 and give a more systematic identification of the scaling. 

Recall that the basic model satisfies 

Xi{t)=Xi{0) + Yb(^J^ K2X2{s)ds^ -Yd (^J^ KiXi {s)ds 

X2{t) = X2(0) + Ya (^J^ KiXlis) ds^ - Yb (^J^ K2X2{s) ds 

-Yf(^J^KeX2{s)Xs{s)ds^, 
X^it) = XsiO) + Y, ^* K^Xiis) ds^ - Ye ^* K^Xsis) ds 
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rt 



With reference to Section 1.2, we consider a general scaling Zl^{t) = A'"~"^Xj(t), 
and replace Kk by XkN^''. Once the Pk are selected, the are determined 
by setting 

for the rate constants Kk given in Section 3 and some appropriate No. 
The normalized system becomes 

rt 



-N-^^Ydl^j\iN^^Z^{s) ds^ , 
Z^{t) = Z^{<d)+N-^Wa{^j\^N^^Z^{s)d.^ 
- iV-°2 Yb (^J^ K2iV°' Z^{s) ds' 



Zi^{t) = Z^{Q)+N-^-'Y^(^j\:iN^^Z^{s)ds 



- N-'^We ^ K^N'^^Z^is) ds 

- N'-^Wf ^* K6Af"^+"»Z2^(s)Zf (s) ds^ . 
Setting Vfc^(t) = Z^{N^t) and replacing by XiN^\ 

V^{t) = Vf(S)) + iV^^in ( r \2N"<+'^m'^Wi^ [s) ds 



- N~'''Yd(^J^ XiN^+'^^N'^Wfis) ds^ , 



iV-"2 Yb ^* A2iV^+^' iV"' ^2^ (s) ds^ 
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We assume that (Fi^(O), ^"^(0), 1/3^(0)) ^ (^1(0), ^2(0), ^3(0)). 

The question is how to determine, in a systematic way, what the expo- 
nents Oi, (3k and 7 should be. There are several conditions that help this 
determination. First, we want the scaling to ensure that Vf^ {t) = 0{\). This 
requirement can be met either by ensuring that the individual terms on the 
right are 0(1) or by ensuring that terms cancel. Second, it is natural to 
assume that the (3k have the same order as the k^, that is, we should have 

(A.6) /?6</32</34</3i</35</33. 

As is clearly reasonable, we assume that (3i = 0. This last assumption is not 
really a restriction, since if (3i / 0, we can add /3i to 7 and substract (3i from 
each of the (3k- 

Finally, comparing the k^, it is also natural to assume that (3^ > (3^ and 
(3q < (32- (We will see that the second of these assumptions is actually implied 
by other considerations.) For the scaling used in Section 3, (3i = (3^ = (3^ = 0, 
(32 = -2/3, P3 = 1 and Pe = -5/3. 

Suppose, as is the case in Section 3, we also require that the scaling makes 
each of the terms in the equation for to be 0(1). In particular, we look 
for a scaling in which the nonlinear behavior is preserved. Then we must 
have 

"2 = 7 + "1, 

a2 = 7 + /?2 + "2, 

"2 = 7 + /?6 + "2 + "S- 

In addition, for V^^ to be 0(1) without being asymptotically negligible, we 
must have 

"2 = 7 + /?2 + "2 = 7 + /34 + ai- 
Similarly, for , we must have 

7 + /Js + ai > 7 + /?5 + as, 

7 + /?3 + "1 > 7 + /^e + "2 + as = "2, 

with equality holding for at least one of the inequalities. Since we are as- 
suming that /Js > /35 > 0, we must have 

7 + /?3 + ai = 7 + /?5 + as 
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and, hence, 03 > ai. 

The above assumptions imply 

Q2 - ai = /34 - /32 = 7 = -/?2 > 0, 

so a2>ai, ^4 = 0, and 

"3 = -1 - P6 = P2- Pg- 

These restrictions leave three cases of interest: ai = 02 > and = 0, 
«! = a2 > and > 0, and a2> ai> 0. 

If ai = 02 > 0, then 02 = Pi = Pa = and as = —Pq. If, in addition, P^ = 0, 
then 03 = ^ + P3 + ai and as — > 00, the system converges to the solution 
of 

Vi{t)=Vi{0) + I\\2V2{S) - \AVi{s))ds, 

Jo 

V2{t) = 1/2(0) + f\xiVi{s) - X2V2{S) - XeV2{s)V3{s)) ds, 
Jo 

Vsit) = VsiO) + f\x3Vi{s) - X^Vsis)) ds. 
Jo 

If = a2 > and /?5 > 0, then 

(A.7) hm ( [\xsVf{s)-X,Vf{s))ds)=0, 

and {Vi^ ,¥2^) converges to a solution of 

Vi{t) = Vi{0) + f\x2V2{s) - X4Vi{s))ds 
Jo 

V2{t) = V2{0) + J^[Ws) - X2V2{S) - ^V2{s)Vi{s)^ ds. 

If a2 > ai, then P2 < and 7 > and, consequently, 7 + /?4 > and 
7 + /Js > 0. It follows that 

(A.8) lim ( f\x4V/^{s) - A2V/(s)) ds] = 

iV->oo\Jo / 

and (A.7) hold. Then, as in the scaling in Section 3, V2^ converges to the 
solution of 

(A.9) V2{t) = V2{0) + /^*((^-A2)F2(.)-^V2(.)^)d.. 

Define Xk = KkN^^'^ for some Nq. Then 

A1A2 / H1H2 \ ,.7 
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and 

A5A4 K5K4 

— 1\q 

Recalling that V^°{t) = NQ°''^X2{N^t), the convergence suggests approxi- 
mating X2it) by V2{t) = Nl^'V2iNQ^t). But if V2 satisfies (A.9), then V2 
satisfies 

(A.IO) y2(t)=V2(0)+ f((^-K2)v2is)-^^V2isf)ds, 

SO the approximation does not depend on the choice of the scaling parameters 
beyond the restrictions identified above and the assumption that a2> ol\. 

The behavior of the V^^ and depends primarily on whether ai > or 
a\ = 0. If ai > 0, then (A. 7) and (A. 8) can be strengthened to 

hm sup ilXsVf'is) - AgKf + |A4^i^(s) - Aal^a'^Cs)]) = 
for each < e < t. 

If ai = 0, then the behavior of V/^ is essentially the same as in Section 3. 
If, in addition, = 0, then the joint behavior of and is essentially 
the same as in Section 3. If ai = and /^s > 0, then 

(A.ll) hm f\X3Vf{s)-X,Vf{s)\ds = 0, 

and for each t > 0, {Vf {t),V3^ {t)) {Vi{t), ^Vi(t)), where Vi{t) is Poisson 
distributed with parameter ^^^^J^^'^ ■ 
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